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INTRODUCTION 


This  is  the  second  in  a  series  of  four  reports  whose  overall  purpose  is  to  describe  an 
adaptive  finite  element  method  (AFEM)  for  solving  systems  of  parabolic  partial  differential 
equations.  In  particular,  AFEM  attempts  to  find  a  numerical  solution  of  an  M-dimensional 
system  of  the  form 

ut(x,t)  +  f(x,t,u,ux)  =  [D(x,t,u)ux(x,t)]x  ,  a<x<b  ,  t> 0,  (la) 

subject  to  the  initial  conditions 

u(x,0)  =  u°(x)  ,  azxzb  (lb> 

and  linear  separated  boundary  conditions 

Al(t)u(a,t)  +  Bl(t)ux(a,t)  =  gl(t)  (lc) 

A\t)u(b,t)  +  BT(t)ux(b,t)  =  g\t )  ’  l> 


The  variables  x  and  t  represent  spatial  and  temporal  coordinates  and  denote  partial 
differentiation  when  they  are  used  as  subscripts;  u,  f,  u°,  gf,  and  g  are  M-vectors;  and  D,  A1,  B,  AT, 
and  R  are  MxM  matrices. 

The  problem  is  assumed  to  be  well-posed  and  parabolic;  thus,  e.g.,  D(x,t,u)  is  positive 
definite.  The  rows  of  B  and  W  are  restricted  to  be  either  entirely  zero  or  a  row  of  the  MxM 
identity  matrix.  When  the  ith  row  of  B  or  &  is  identically  zero,  then  A~  or  Ad  cannot  be  zero, 
respectively,  and  the  boundary  condition  is  a  Dirichlet  (essential)  condition.  Otherwise,  the 
boundary  condition  is  a  Neumann  or  Robbins  (natural)  condition.  The  ultimate  goal  of  AFEM 
is  to  determine  an  approximate  solution  to  Eq.  (1)  to  within  a  user  prescribed  error  tolerance. 

The  adaptive  strategies  utilized  by  AFEM  are  (1)  error  estimation  coupled  with  (2)  local 
mesh  refinement  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  1),  Babuska  and  Dorr  (ref  2),  Babuska, 
Zienkiewicz,  Gago,  and  Oliveira  (ref  3),  Bank  and  Weiser  (ref  4),  Berger  and  Oliger  (ref  5), 
Bieterman  and  Babuska  (refs  6,7),  Moore  and  Flaherty  (ref  8),  Shephard  (ref  9),  Strouboulis  and 
Oden  (ref  10),  Zienkiewicz  and  Zhu  (ref  11)),  and  (3)  mesh  movement  (cf.,  e.g.,  Adjerid  and 
Flaherty  (ref  1),  Amey  and  Flaherty  (ref  12),  Bell  and  Shubin  (ref  13),  Davis  and  Flaherty  (ref 
14),  Dorfi  and  Drury  (ref  15),  Dwyer  (ref  16),  Ewing,  Russell,  and  Wheeler  (ref  17),  Hyman  (ref 
18),  Kansa,  Morgan,  and  Morris  (ref  19),  Miller  and  Miller  (refs  20,21),  Petzold  (ref  22),  Rai  and 
Anderson  (ref  23),  Russell  and  Ren  (ref  24),  Saltzman  and  Brackbill  (ref  25),  Smooke  and 
Koszykowski  (ref  26),  Thompson  (ref  27),  Verwer,  Blom,  Furzeland,  and  Zegeling  (ref  28),  and 
White  (ref  29)). 
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The  purpose  of  this  report  is  to  describe  the  error  estimating  procedures  employed  by 
AFEM.  Detailed  summaries  of  how  AFEM  implements  its  other  adaptive  strategies  are  found  in 
separate  reports  entitled:  Adaptive  Finite  Element  Method  III:  Mesh  Refinement  (ref  30)  and 
Adaptive  Finite  Element  Method  IV:  Mesh  Movement  (ref  31).  Furthermore,  the  report, 
Adaptive  Finite  Element  Method  I:  Solution  Algorithm  and  Computational  Examples  (ref  32), 
describes  how  all  the  adaptive  algorithms  are  implemented  in  unison  and  contains  results 
demonstrating  the  utility  of  AFEM  as  a  computational  tool. 

The  error  estimation  performed  by  AFEM  is  based  on  the  work  of  Adjerid  and  Flaherty 
(ref  1).  Adjerid  and  Flaherty  developed  an  a  posteriori  estimate  to  the  spatial  discretization 
error  of  a  finite  element  method  of  lines  solution  for  a  vector  system  of  parabolic  partial 
differential  equations.  They  discretized  the  system  in  space  using  Gaierkin’s  method  with 
piecewise  polynomial  finite  element  approximations  of  an  arbitrary  order  p.  The  error  estimate 
was  calculated  using  Gaierkin’s  method  with  piecewise  polynomial  functions  of  order  p  +  1.  A 
nodal  superconvergence  property  of  the  finite  element  method  was  used  to  neglect  errors  at 
nodes,  and  thus,  improve  computational  efficiency.  Ordinary  differential  equations  (ODEs)  for 
the  finite  element  solution  and  error  estimation  were  then  integrated  in  time  using  the  backward 
difference  code  DASSL  (cf.,  Petzold  (ref  33)). 

Adjerid  and  Flaherty  (ref  1)  assumed  that  the  temporal  discretization  error  associated 
with  DASSL  was  negligible  compared  to  the  spatial  error.  Thus,  their  estimate  of  the  spatial 
discretization  error  could  be  regarded  as  an  estimate  of  the  total  error.  They  used  their  error 
estimate  to  control  mesh  moving  and  local  mesh  refinement  procedures  that  simultaneously 
attempted  to  equidistribute  the  error  estimate  and  satisfy  a  prescribed  global  error  tolerance. 
Similar  mesh  refinement  strategies  have  been  used  by  Bieterman  and  Babuska  (refs  6,7). 

Our  goal  is  to  develop  techniques  that  simultaneously  estimate  temporal  and  spatial 
discretization  errors.  With  such  estimates,  mesh  refinement  and/or  moving  decisions  can  be 
made  to  reduce  the  largest  component  of  the  error  with  the  least  amount  of  work.  For  example, 
if  the  temporal  error  is  the  dominant  component  of  the  total  error,  then  one  need  only  adjust  the 
time  step  in  order  to  improve  accuracy.  In  this  way,  one  avoids  needlessly  increasing  the  spatial 
discretization  which  would  increase  the  computational  complexity  unnecessarily.  Local  and 
global  estimates  of  the  discretization  error  have  been  successfully  used  to  control  refinement 
algorithms  that  attempt  to  solve  partial  differential  equations  to  prescribed  levels  of  accuracy  (cf., 
e.g.,  Babuska,  Zienkiewicz,  Gago,  and  Oliveira  (ref  3)  and  Flaherty,  Paslow,  Shephard,  and 
Vasilakis  (ref  34)  for  a  sampling). 

As  in  Adjerid  and  Flaherty  (ref  1),  Eq.  (la)  is  discretized  in  space  using  Gaierkin’s 
method  with  piecewise  linear  finite  elements.  Temporal  discretization,  however,  is  performed  by 
the  backward  Euler  method  as  opposed  to  using  an  ODE  code  (cf.,  Coyle  and  Flaherty  (ref  32)). 
A  second  solution  is  calculated  using  trapezoidal  rule  integration  in  time  and  the  difference 
between  the  two  solutions  is  used  to  furnish  an  estimate  of  the  temporal  discretization  error.  A 
third  solution  is  obtained  using  quadratic  finite  elements  and  the  trapezoidal  rule  in  time.  This 
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solution  is  higher  order  in  space  and  time  than  the  original  piecewise  linear  finite  element- 
backward  Euler  solution.  Hence,  it  can  be  used  to  provide  an  estimate  of  the  total  discretization 
error  of  the  piecewise  linear  finite-element  backward  Euler  solution.  Furthermore,  the 
difference  between  the  piecewise  linear  and  quadratic  solutions  calculated  by  the  trapezoidal  rule 
furnishes  an  estimate  of  the  spatial  discretization  error  (cf.,  Moore  and  Flaherty  (ref  8)  or  Coyle 
and  Flaherty  (ref  35)). 

At  first  sight,  the  above  procedure  seems  expensive;  however,  nodal  superconvergence 
significantly  reduces  computational  complexity.  In  the  present  context,  superconvergence  implies 
that  finite  element  solutions  converge  at  a  faster  rate  at  mesh  point  locations  than  elsewhere  in 
the  problem  domain  (cf.,  Adjerid  and  Flaherty  (ref  1)).  Hence,  the  error  at  the  nodes  may  be 
neglected  relative  to  the  error  in  the  interior  of  the  elements  when  N,  the  number  of  mesh 
points,  is  sufficiently  large.  Nodal  superconvergence  has  been  used  by  several  investigators  as  a 
means  of  constructing  a  posteriori  error  estimates  in  finite  element  approximations  (cf.,  Adjerid 
and  Flaherty  (ref  1),  Bieterman  and  Babuska  (refs  6,7),  and  Coyle  and  Flaherty  (ref  12)).  Defect 
correction  methods  can  also  be  used  to  reduce  costs  associated  with  the  temporal  integration  (cf., 
Dahlquist,  Bjork,  and  Anderson  (ref  36)). 

The  piecewise  linear  and  quadratic  finite  element  procedures  and  the  temporal 
integration  schemes  are  outlined  in  the  Numerical  Discretization  section  (cf.,  Coyle  and  Flaherty 
(ref  32)  for  a  more  complete  description).  Derivations  of  the  various  error  estimates  (total, 
spatial,  and  temporal)  are  presented  in  the  Error  Decomposition  and  Estimates  section.  Then  in 
Convergence  Examples,  examples  that  indicate  the  convergence  of  the  error  estimates  to  the  true 
error  and  its  components  are  described.  Finally,  in  the  last  section,  a  summary  of  this  report  is 
presented. 


NUMERICAL  DISCRETIZATION 

A  weak  form  of  the  problem  is  constructed  by  multiplying  Eq.  (la)  by  a  test  function 
v(x,t)  G  Hg ,  integrating  the  result  with  respect  to  x  from  a  to  b,  and  integrating  the  diffusive 
term  by  parts  to  obtain 

(v,ut)  +  (vj)  +  A(v,u)  =  vTDux\ba  ,  t> 0  ,  for  all  v  e  h£  .  ( 


The  inner  product  (v,u)  and  strain  energy  A(v, u)  are  defined  as 

(v,u)  =  fbvrudx,  A(v,u )  =  jbvTxDux  dx  .  (2b, c) 

Functions  v  belonging  to  H1  are  required  to  have  finite  values  of  (v,v)  and  (vx,vj.  Functions  in 
Hi  are  in  H1  and  vanish  at  x  =  a  and/or  b  if  an  essential  boundary  condition  is  applied  there. 
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Any  weak  solution  u  G  H2E  of  Eq.  (2a)  must  also  satisfy  any  essential  boundary  conditions  at  x  = 
a 


or  at  x  =  b 


M 

Y,A\ft)ufa,t) 


j*i 


*  a'« 


M 


*  A^t) 


1*1 

;=i 


(2d) 


(2e) 


when  the  ih  row  of  B1  and/or  B  is  zero,  respectively.  Natural  boundary  conditions  replace  the  ith 
component  of  ux  at  x  =  a  or  b  in  Eq.  (2a)  when  prescribed. 

Initial  conditions  for  Eq.  (2a)  are  obtained  by  L2  projection,  i.e., 

(v,u)  =  (v,«°)  ,  t  =  0  ,  for  all  v  e  Hq  .  ^ 


A  discrete  version  of  the  weak  system  Eq.  (2)  is  constructed  by  using  finite  element- 
Galerkin  procedures  in  space  and  finite  difference  techniques  in  time  on  a  fully  adaptive  mesh 
(one  that  is  both  refined  and  moved  as  time  progresses). 

Spatial  Discretization 

To  discretize  Eq.  (2a)  in  space,  introduce  a  time-dependent  partition 

H^f)  =  {  a  =  x0<xl(t)<x2(t)<...<xN  =  b  }  (3 

of  (a,b)  into  N  subintervals  (xu(t),  xft)),  i=l,2,...,N  and  approximate  u  and  v  by  piecewise 
polynomial  functions  U  and  V,  respectively,  with  respect  to  this  partition.  Thus,  the  spatially- 
discrete  form  of  Eq.  (2a)  consists  of  finding  U  G  C  H2E  such  that 
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(4a) 


(v,ut)  +  (vj)  +  i4(y,i/)  =  vtdux  i*  ,  r>o , 

/or  a//  V  €  Sq  Hq  , 


(V,U)  =  (V,u°)  ,  t= 0  ,  for  all  V  e  S0N  c  flj  .  (4b) 


The  spaces  S£  and  5^  will  consist  of  either  piecewise  linear  or  piecewise  quadratic  polynomial 
functions.  The  spaces  of  piecewise  linear  polynomials  are  denoted  S^’1  and  S^1  and  a  basis  is 
easily  constructed  in  terms  of  the  familiar  "hat"  functions 


ft) 

xft)  -x^ft) 
Xi+1  (t)-x 

Xi+1  ft)  -xft) 

0 


Xj-xd)  <  X  z  xft) 
xft)  <S  X  z  Xi+1(t)  . 


otherwise 


(5) 


The  piecewise  linear  finite  element  solution  U2  G  S^1  is  written  in  the  form 

N 

i= 0 


and  determined  by  solving  the  ordinary  differential  system 


(VvU1()  +  (Vvf(%t,UvUlx))  +  A(VvUx)  =  V\DUlx  t  , 


f>0  ,  for  all  Vj  e  S0 


N,1 


(7a) 


(V^)  =  (Vvu°)  ,  r  =  0  ,  for  allVl  e  S0W  , 


(7b) 


where  the  piecewise  linear  test  functions  V1  G  SN0J  have  a  form  similar  to  Eq.  (6). 
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Piecewise  quadratic  approximations  V2  G  SN0’2  are  constructed  by  adding  a  "hierarchical" 
correction  E2(x,t)  to  U},  i.e., 


U2(x,t )  =  Ux(x,t)  +  E2(x,t)  , 


where 


N 


E2{xj)  =  Y,di-vP)fi-v2(xt)  • 


i=l 


(8a) 


(8b) 


The  basis  i=l,2,...,N,  for  the  quadratic  correction  has  the  form 

2[x-xi_1(f)][x-xl(f)] 


=  \ 


0 


,  XM(f)  <  X  <  x.(t) 


otherwise 


Piecewise  quadratic  solutions  are  determined  by  solving 


(V2,U2{)  +  (V2,f(\t,U2,U2))  +  A(VvU2)  =  iJ  , 


t> 0  ,/ora//  y2  6  S™  , 


(9) 


(10a) 


(V2,U2)  =  (V2,u°)  ,t= 0,  for  all  V2  e  S™  , 


(10b) 


where  V2  has  a  form  similar  to  Eq.  (8). 

Temporal  Discretization 

Discretization  in  time  is  performed  by  integrating,  for  example,  Eq.  (4a)  over  the  time 
step  from  f'1  to  f  to  obtain 
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(11a) 


N 


for  all  V  €  Sq  , 


(V,U)  =  (V,u°)  ,  t=0  ,  for  all  V  e  S0N  .  (H*>) 


The  integration  in  Eq.  (11)  will  be  over  an  irregular  region  due  to  the  mesh  motion  with  the  test 
function  V  having  an  undesirable  time  dependency. 


In  order  to  overcome  this  difficulty,  introduce  a  linear  transformation 
*  =  Xm1  +  +  VfcAx^l+O  +  V2(Ax,n-Axrl)v(l^), 


(12a) 


t  =  f"_1  +  AtnT  (12b) 

where 

Ajc"  =  Jt,n  -  ,  Af"  =  tn  -  tn~x  (12c»d) 

and 

x-  =  xft")  ,  i  =  0,1,...^  .  (12e) 


This  transformation  maps  the  rectangle  {(£,t)|-1  <f<l,  0<t<1}  onto  the  trapezoidal 
space-time  element  whose  vertices  are 

.  tfV')  ,  (sitjf)  and  (*,",(")  • 


The  basis  elements  $iA(x,t)  and  fa (x,t),  the  only  nonzero  ones  on 
{(x,t)  | \xhl(t)  <  x  <  xft) ,  f1  <  t  <  f},  are  transformed  to  functions  with  no  r  dependency;  thus, 
4>i-i(x,t)  and  fafct)  become,  respectively, 
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*_!«)  =  V2(l-0  ,  and 


(13a) 


=  V2(i+o  ,  -l  *5  *  l  .  d3b) 

Define 

F.  =  J*1  (yrl/,  +  yr/  +  vjDiy  dx  dt  (14) 

and  write  Eq.  (11)  as 

V  F-  =  f1  yrWj*  AfVx  ,  for  all  V  e  S0N  ,  d5a) 

i= i  Jo 

(yt/)  =  (y,K°)  ,  t= 0  ,/or  all  V  €  Sq  .  (15b) 


by  performing  the  change  of  variables  from  t  to  r  (cf.,  Eq.  (12b)). 

Transforming  Eq.  (14)  from  the  (x-t)- plane  to  the  (£,r)-plane  (cf.,  Eq.  (12))  yields 


■IX 


Vr(Uxl)x  -  VT(Ux\tx  +  VTJx^tx  +  vlDUl  — 

H 


dldx 


(16a) 


where 


x  = 


(16b) 


Equation  (16a)  can  be  simplified  further  by  integrating  by  parts  to  obtain 

F.  =  G,(l)  -  G.(0)  +  A tn[1  If t)  dx  (17a) 

11  1  Jo 

where 


8 


G,(t)  -  /'  rtfc,  <«  , 


(17b) 


/,M  =  /'  [-Vr(£tt)f  +  Vr/t{  + 


(17c) 


Substituting  Eq.  (17a)  into  Eq.  (15a)  then  yields 

AT 


£  G.(l)  -  G/0)  +  At"  f 1  //T)*]  =  Arf1  VTDU  \ba  dx  , 

j_j  **0  «<0 


for  all  V  €  Sq  , 


(18a) 


(V,U)  =  (V,«°)  ,  t=0  ,  /or  o//  V  e  Sq  . 


(18b) 


All  that  remains  is  to  approximate  the  time  integrals  in  Eq.  (18)  using  a  quadrature  rule. 
This  is  done  by  using  a  weighted  two-step  method,  which  for  Eq.  (18)  has  the  form 


AT 


£  [G.(l)  +  G.(0)  +  At "07.(1)  +  At "(1-0)7.(O)]  =  At*dVTDUx\ba\x=l 
+  At”(l -ff)VTDUx\ba |t=0  Jorall  VeS0N,de  [0,1] , 


(19a) 


(V,U)  =  (y,u°)  ,t=0,  for  all  V  e  S0N  . 


(19b) 


Only  two  choices  of  6  are  utilized:  either  0  =  1,  which  yields  the  backward  Euler  method,  or  0  = 
V2,  which  yields  the  trapezoidal  rule. 


ERROR  DECOMPOSITION  AND  ESTIMATES 

Strang  and  Fix  (ref  37)  demonstrate  how  the  total  discretization  error  of  a  numerical 
method  for  solving  parabolic  partial  differential  equations,  which  couples  finite  elements  in  space 
with  finite  differences  in  time,  can  be  decomposed  into  its  spatial  and  temporal  components. 
They  do  this  by  defining  an  intermediate  solution,  U,  that  satisfies  the  finite  element 
discretization  (cf.,  Eq.  (4a))  but  is  integrated  exactly  in  time.  Continuing  to  let  u  denote  the 
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exact  solution  of  the  Galerkin  problem  (cf.,  Eq.  (2)),  let  Ue  denote  the  numerical  approximation 
which  satisfies  both  the  finite  element  equations,  and  the  finite  difference  equations  in  time  (cf., 
Eq.  (19)).  Then  the  total  discretization  error,  e,  where 

e  =  u  -  UQ  (20a) 


can  be  bounded  as  follows: 

Ik II  =  \u-U+U  -  UJ  *  flu  -  U\\  +  \U  -  tf0||  .  (20b) 


Strang  and  Fix  (ref  37)  show  how  the  first  term  of  the  right-hand  side  of  Eq.  (20b), 

||u  - 1/||,  is  strictly  an  error  due  to  the  finite  element  approximation  process,  and,  as  such, 
dependent  only  on  the  spatial  discretization.  Thus,  |u  -  l/fl  is  the  spatial  component  of  the  total 
discretization  error.  Similarly,  they  show  that  ||t/  -  Ue\\  is  dependent  on  the  temporal 
discretization  and  hence  represents  the  temporal  component  of  the  total  discretization  error. 

Our  goal  is  to  estimate  the  discretization  error  per  time  step  in  solutions  of  Eq.  (19) 
obtained  by  using  piecewise  linear  finite  element  approximations  in  space  and  the  backward 
Euler  method  in  time.  It  seems  most  appropriate  to  gage  errors  in  the  H1  norm 

Nil,  =  K<vO  +  Mf* ;  (21) 


however,  other  metrics  may  also  be  used.  An  error  estimate  that  is  global  in  space  and  local  in 
time  may  at  first  seem  unusual,  but  it  is  commonly  used  when  spatial  finite  element 
approximations  are  combined  with  temporal  finite  difference  methods  (cf.,  Thom<5e  (ref  38)). 

Let  the  piecewise  linear  finite  element  solutions  at  time  f1  obtained  by  using  backward 
Euler  ( 0=1  in  Eq.  (19))  and  trapezoidal  rule  (9=V. 2  in  Eq.  (19))  temporal  integration  be  denoted 
by  U"  and  U?/2,  respectively.  Likewise,  let  U?/2  denote  the  piecewise  quadratic  finite  element 
solution  of  Eq.  (19)  at  f  found  by  using  the  trapezoidal  rule  integration  in  time. 

It  is  known  (cf.,  Strang  and  Fix  (ref  37))  that 

II -  ^(Ollj  =  0((Af*)2)  +  OtAT1)  .  (22> 


Since 


WV*>  -  -  two3)  +  0(N-2) , 


(») 


one  should  be  able  to  use  the  difference  between  U?/2  and  U\  to  estimate  the  error  in  U] ;  thus, 
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(24) 


ii**  ‘•'ini  11"%  "i hi  •  ii"  "%ni 

s  |0£  -  I^llj  +  0((A»”)3)  +  W2)  . 

The  main  problem  in  using  \U?A  -  V}  ^  as  an  a  posteriori  estimate  of  (| u  -  V}  I,  is  the 
computational  effort  required  to  obtain  V?A.  This  cost  can  be  reduced  considerably  by  using  the 
superconvergence  property  of  the  finite  element  method  for  one-dimensional  parabolic  systems. 
Nodal  superconvergence  enables  us  to  approximate  U?A  as 

K  -  K  *  K  (25) 


where  U?A  is  obtained  by  solving  Eq.  (19)  using  trapezoidal  rule  integration  and  E?A  is  obtained  by 
solving  Eq.  (19)  by  trapezoidal  rule  integration  \yjth  Z7  replaced  by  Eq.  (8b).  Furthermore,  it  is 
only  necessary  to  test  Eq.  (19)  against  functions  V  E  SN0’2,  where  SN0-2  is  a  space  of  quadratic 
polynomials  that  vanish  on  IlN(t). 

To  summarize,  the  procedure  for  obtaining  the  finite  element  solution  lf\  and  its  error 
estimate  U?A  +  E?A  -  U}  for  the  time  step  [f1,?]  consists  of: 

1.  Determining  U}  as  the  solution  of 

E[G,(1)  -  *  A<*/((1)]  -  A , 

i= 1 

for  all  V  £  Sq’1 

where 

G{1)  =  f  \  VTU?X*  d\  ,  Gf 0)  =  vTif;xxn{1  dt , 


(26a) 


(26b, c) 


w  -  /: 


-VT{Ulx\  +  VTrn  +  V^DnUnl 


l  n 


dl  , 


(26d) 


and 


(V,l$  =  (V,u°)  ,t= 0  ,  for  all  V  e  S*1  . 


(26e) 


2.  Determining  V?/2  as  the  solution  of 

n  r 

E  Gfx>  -  g,<°)  +  Ko) +  W  = 

i=l  2 

-^f VTDnU?A  +  VTDn-ilfl'x  ]b  ,for  all  V  e  S*’1 

2  *  A  a 

where 

G,d)  *  VV^c’  df,  ,  (3,(0)  =  VT^-'x’-1  d\  , 

/,(!)  =  /’  -V(l4i")s  +  VTJ”x”  +  d5 , 

. 

/(°)  =  /  ^  -yr(i^‘1in-1)€  +  y^r1  +  , 

and 

(y,t^)  =  (y,«°)  ,  *=0  ,  /or  atf  y  e  So""’1  . 

3.  Determining  E"2  as  the  solution  of 

E  d.w  -  ^(°> +  «•  t,(o>]  -  o , 

i= 1  L  ** 

for  all  V  g  S** 

where 

G/D  -  /'  Vr£.X«  ,  G,(0)  =  pJT^r'di, 

and 


(27a) 

(27b, c) 

(27d) 

(27e) 

(27f) 

(28a) 

(28b,c) 
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(28d) 


(28e) 

(28f) 


Temporal  error  estimation  is  local;  thus,  we  use  VJ1  as  an  initial  condition  for  the 
trapezoidal  rule  integrations  in  Eqs.  (27)  and  (28).  Nodal  superconvergence  and  the  hierarchical 
formulation  has  uncoupled  the  piecewise  linear  and  quadratic  components  of  U?A.  The  spatial 
error  estimate  E?A  on  the  subinterval  (x^,  x;)  is  furthermore  uncoupled  from  the  error  on  other 
subintervals  and  this  significantly  reduces  the  computational  complexity  associated  with  solving 
Eq.  (28).  The  solution  of  Eq.  (27),  noted  in  Step  2,  is  necessary  in  order  to  increase  the 
temporal  accuracy  of  the  solution  because  superconvergence  only  increases  the  order  of  accuracy 
in  space.  Some  computational  savings  can  generally  be  obtained,  especially  for  nonlinear 
problems,  by  calculating  V?A  as  a  defect  correction  to  the  backward  Euler  solution  lFr 

As  described  above, 


K 


-  «?!, 


(29) 


furnishes  an  estimate  to  the  error  ||u  -  U\  lx  of  the  backward  Euler  solution.  Equation  (29) 
suggests  the  inequality 


en  s  || -  iflr  +  iKlIi  . 


(30) 


The  term  \U?A  -  U] |t  is  the  difference  between  two  piecewise  linear  solutions  computed  with 
temporal  integration  schemes  of  different  orders  and  can  be  regarded  as  a  measure  of  the 
temporal  discretization  error.  In  a  similar  manner,  ||£"J1  can  be  regarded  as  a  measure  of  the 
spatial  discretization  error.  Indeed,  when  the  finite  element  system,  Eq.  (4),  is  integrated  exactly 
in  time,  Adjerid  and  Flaherty  (ref  1)  proved  that  |E||j  converges  to  the  exact  spatial 
discretization  error  Hu  -  U Ih  as  N  ->  «>  for  linear  parabolic  problems. 
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CONVERGENCE  EXAMPLES 


Example  1 


Consider  the  linear  heat  conduction  problem 

,  0  <  x  <  1  ,  f>0. 

(31a) 

u(x, 0)  =  simtx:  ,  0  ^  x  ^  1  , 

(31b) 

u(0,t)  =  u(l,t)  =  0  ,  t  >  0  . 

(31c, d) 

The  exact  solution  of  this  simple  problem  is 

u(x,t)  =  e sin  ttx  . 

(32) 

We  solved  Eq.  (31)  on  a  uniform  mesh  with  N  finite  elements  for  one  time  step  At  using 
the  methods  described  above  and  several  choices  of  N  and  At.  The  effectivity  index 

e  :  = - - - —  (33) 

W-.A 0  -  I'l  l, 

(cf.,  Babuska,  Miller,  and  Vogelius  (ref  39)),  is  used  as  a  means  of  gaging  the  accuracy  of  the 
error  estimate  e1.  Ideally,  we  would  like  ©  not  to  differ  appreciably  from  unity  and  to  approach 
unity  as  jV  and  At  -»  0. 

We  present  a  summary  of  results  for  a  sequence  of  calculations  performed  with  N  —  2™ 
and  At  =  1.024x2m,  m  =  1, in  Figure  1.  These  results  strongly  suggest  that  ©  -*  1  as  m  -» 
00. 
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Total  Effectivity  vs.  m 

Number  of  elements  =  2W 
Time  step  =  1.024  x  Tm 


-i - 1 - 1 _ i _ i _ i _ i _ i _ i 

23456789  1( 


m 


1.  Total  eff< 


discretization  for  Example  1. 


Example  2 


Consider  the  linear  heat  conduction  problem 

ut  +  u  =  —  ,  0  <  x  <  1  ,  t  >  0  , 

n2 


(34a) 


u(x, 0)  =  1  ,  0  ^  x  <,  1  , 


(34b) 


u(0,t )  =  =  e  ‘  ,  t  >  0  . 

The  exact  solution  of  this  simple  problem  is 

u(x,t)  =  e  ‘  . 


(34c, d) 


(35) 


We  solved  Eq.  (34)  on  a  uniform  mesh  with  eight  finite  elements  for  one  time  step  At 
using  the  methods  described  above  and  several  choices  of  At.  The  temporal  effectivity  index 


0.  := 


\K  -  ulh 

ll«(‘,Af)  -  UX 


(36) 


is  used  as  a  means  of  gaging  the  accuracy  of  the  error  estimate  ]| U}/2  -  U\  ||r  Again,  we  would 
like  @t  not  to  differ  appreciably  from  unity  and  to  approach  unity  as  At  -»  0,  since  there  is  no 
spatial  component  to  the  total  discretization  error. 

We  present  a  summary  of  results  for  a  sequence  of  calculations  performed  with  At  = 
1.024x2m,  m  =  in  Figure  2.  These  results  strongly  suggest  that  ©t  -»  1  as  m  -*  °o. 
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Example  3 


Consider  the  linear  heat  conduction  problem 


1 

ut  -  u  =  -r  “«  » 

0<x<l,t>0, 

(37a) 

7t 

u(x, 0)  =  sin  tcjc 

,  0  ^  x  <>  1  , 

(37b) 

u(0,t)  =  u(\,t) 

=  0  ,  t  >  0  . 

(37c, d) 

The  exact  solution  of  this  simple  problem  is 

u(x,t)  = 

sin  tix  . 

(38) 

We  solved  Eq.  (37)  on  a  uniform  mesh  with  N  finite  elements  for  one  time  step  A t  = 
0.001  using  the  methods  described  above  and  several  choices  of  N.  The  spatial  effectivity  index 

e  :=  ...  (39) 

|«(*,Af)  -  UX 

is  used  as  a  means  of  gaging  the  accuracy  of  the  error  estimate  ||2?4li*  Again,  ©s  should  not 
differ  appreciable  from  unity  since  for  such  a  small  At,  there  is,  effectively,  no  temporal 
component  to  the  total  discretization  error. 

We  present  a  summary  of  results  for  a  sequence  of  calculations  performed  with  N  =  2", 
m  =  in  Figure  3. 
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SUMMARY 


Methods  for  calculating  a  posteriori  estimates  of  the  total,  spatial,  and  temporal 
discretization  errors  when  a  vector  system  of  parabolic  partial  differential  equations  is  solved 
using  piecewise  linear  finite  elements  in  space  and  the  backward  Euler  method  in  time  was 
presented.  First,  the  division  of  the  total  discretization  error  into  its  spatial  and  temporal 
components  was  shown  theoretically.  This  was  followed  by  a  method  to  approximate  these  errors 
numerically.  Then  it  was  shown  how  to  obtain  these  error  estimates  by  using  higher-order 
methods,  with  nodal  superconvergence,  in  order  to  improve  computational  efficiency.  Finally,  a 
comparison  of  the  exact  and  estimated  errors  was  presented  in  Examples  1,  2,  and  3  and  in 
Figures  1,  2,  and  3. 

Comparison  of  the  exact  and  estimated  errors,  presented  in  Examples  1,  2,  and  3,  give  us 
some  confidence  in  the  accuracy  of  our  error  estimates.  As  indicated  by  Figures  1,  2,  and  3,  the 
three  estimates  all  converge  to  the  true  errors  as  the  appropriate  discretization  levels  are 
increased.  Even  for  coarse  levels,  results  indicate  that  the  estimates  do  a  reasonable  job  of 
approximating  the  exact  error  (cf..  Figures  1,  2,  and  3).  Thus  not  only  does  one  get  an  indication 
when  the  total  error  of  the  method  is  too  large,  but  also  the  dominant  component  of  the  error. 
With  this  knowledge,  one  can  adjust  the  appropriate  discretization  level  accordingly.  Even  in 
situations  when  the  estimates  are  not  accurate,  one  gets  an  indication  that  refinement  is 
necessary  in  a  particular  region. 
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